function F = fred(k,v,onevar)
% FRED  Fredholm integral operator.
% 
% F = FRED(K,V) computes the Fredholm integral with kernel K:
%    
%      (F*v)(x) = int( K(x,y)*v(y), y=a..b ),
% 
% where [a b] = domain(V). The kernel function K(x,y) should be smooth for 
% best results.
%
% K must be defined as a function of two inputs X and Y. These may be
% scalar and vector, or they may be matrices defined by NDGRID to represent
% a tensor product of points in DxD. 
%
% FRED(K,V,'onevar') will avoid calling K with tensor product matrices X 
% and Y. Instead, the kernel function K should interpret a call K(x) as 
% a vector x defining the tensor product grid. This format allows a 
% separable or sparse representation for increased efficiency in
% some cases.
%
% See also domain/fred, chebfun/volt, chebop.

% Copyright 2011 by The University of Oxford and The Chebfun Developers. 
% See http://www.maths.ox.ac.uk/chebfun/ for Chebfun information.

    % Require two inputs.
    if nargin == 1
        error('CHEBFUN:FRED:nargin','Not enough input arguments.');
    end

    % Inputs in correct order. let this slide...
    if isa(k,'chebfun'),  tmp = v; v = k; k = tmp; end

    % Default onevar to false
    if nargin==2, onevar=false; end     

    % Loop for quasimatrix support
    F = chebfun;
    for j = 1:numel(v)
        F(j) = fred_col(k,v(j),onevar);
    end

end

function F = fred_col(k,v,onevar)
    % At each x, do an adaptive quadrature.
    % Result can be resolved relative to norm(u). (For instance, if the
    % kernel is nearly zero by cancellation on the interval, don't try to
    % resolve it relative to its own scale.) 
    v = set(v,'funreturn',0);
    nrmf = norm(v);
    d = domain(v);
    opt = {'resampling',false,'splitting',true,'exps',[0 0],'scale',nrmf};
    int = @(x) sum(chebfun(@(y) feval(v,y).*k(x,y),d,opt{:}));
    F = chebfun( int, d,'sampletest',false,'resampling',false,'exps',[0 0],'vectorize','scale',nrmf);
    F.jacobian = anon(['[Jvu nonConst] = diff(v,u,''linop'');',...
                       'der = fred(k,d,onevar)*Jvu;'],...
                        {'k','v','d','onevar'},{k,v,d,onevar},1,'fred');
    F.ID = newIDnum;
end % fred_col